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ABSTRACT 

The  prediction  of  unsteady  flow  field  in  turbine  blades  as  well  as  in  the  turbomachinery  stages  is 
now  an  affordable  item,  and  is  required  by  the  reduced  margin  for  increasing  efficiency,  stability 
and  life  of  propulsion  components.  The  numerical  tools  are  now  capable  to  run  within  reasonable 
time  3D  unsteady  calculation  for  full  stage,  and  the  new  techniques  on  the  computation  and  parallel 
computer  allow  the  improvements  of  results  in  terms  of  cost  and  accuracy.  Despite  these  advantages 
many  questions  remain  open  and  the  physical  modelling  joint  with  the  numerical  improvements  is 
still  a challenge  if  it  has  to  produce  usable  results,  compared  with  the  experiments.  On  the  other  side 
the  huge  amount  of  data  extracted  from  experiments  require  care  and  skinless  to  become  useful 
tools  for  design.  The  two  activities  interact  and  support  each  other  in  the  attempt  to  improve  design 
quality.  Aim  of  this  paper  is  the  report  on  some  experience  and  the  attempt  to  give  some  answer  on 
that  challenge,  presenting  results  of  a recent  activity  on  modelling  side  compared  with  experiments 
as  well.  A full-3D  unstructured  solver  based  on  an  upwind  TVD  finite  volume  scheme  has  been 
developed  and  applied  to  the  simulation  of  an  unsteady  turbine  stage.  The  development  of  the 
numerical  strategy  is  discussed  with  particular  concern  on  the  validation  of  the  unsteady  model 
through  a comparison  against  experiments,  NISRE  approach  and  a 3D  steady  stage  computation. 
The  present  work  considers  the  application  of  the  fully  unstructured  hybrid  solver  for  internal 
viscous  flows,  as  well.  The  multiblock  version  of  the  solver  developed  for  turbine  is  considered, 
because  of  the  highly  improved  performance  as  compared  to  the  single  domain  version  of  the  code. 
Moreover,  the  high  numerical  costs  involved  in  3D  unsteady  computations  required  the 
development  of  a new  parallel  single  program  multiple-data  version  of  the  numerical  solver.  The 
results  compare  favourably  with  a set  of  time  averaged  and  unsteady  experimental  data  available  for 
the  turbine  stage  under  investigation,  which  is  representative  of  a wide  class  of  aero-engines. 

Nomenclature 


fg’h 

Convective  fluxes 

Greek  symbols 

f 

Frequency 

P 

Density 

J 

Implicit  residual  Jacobian 

T 

numerical  time 

Q 

Conserved  variables  vector 

(0 

Rotational  speed 

p 

Static  pressure 

Subsripts 

PO,TO 

Total  pressure  and  temperature 

1,2,3 

NGV  inlet,  interstage  and  rotor 

R 

Residual  vector 

outlet  planes 

S 

Source  term  residual 

Acronymom 

t 

Physical  time 

GMRES 

Generalised  Minimal  Residual 

V 

Element  volume 

Method 

s 

Curvilinear  abscissa 

ILU 

Incomplete  LU  factorisation 

W 

Mass  flow  rate 

NISRE 

Not  lsoentropic  Simple  Radial 

x,y,z 

Space  coordinate 

Equilibrium 

W 

Mass  flow  rate 

PE 

Processing  Element  of  a Parallel 
Architecture 

Paper  presented  at  the  RTO  AVT  Symposium  on  “Reduction  of  Military  Vehicle 
Acquisition  Time  and  Cost  through  Advanced  Modelling  and  Virtual  Simulation  ”, 
held  in  Paris,  France,  22-25  April  2002,  and  published  in  RTO-MP-089. 
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1.  INTRODUCTION 

Computational  Fluid  Dynamics  (CFD)  has  become  an  effective  tool  in  the  analysis  of 
complex  flows  and  the  design  of  more  efficient  machinery  components,  thanks  to  its  versatility  in 
the  investigation  of  different  working  conditions  and  to  the  capability  in  analysing  overall  and 
detailed  information  about  the  flow.  The  simulation  of  a transonic  turbine  stage  requires  both 
efficiency  and  accuracy  in  order  to  forecast  the  stage  performances  and  the  realistic  representation 
of  the  unsteady  stator/rotor  interaction.  In  this  regard  many  applications  of  structured  codes  for  3D 
turbine  investigations  have  been  reported  in  literature. 

A basic  assessment  considering  the  research  activity  and  the  main  impact  of  unsteadiness 
phenomena  is  given  by  Sharma  et  al.,  (1992).  Among  all  the  approaches  proposed  some  authors  use 
implicit  ADI  factorisation  schemes  for  the  simulation  of  quasi-3D  stator-rotor  interaction,  while 
others  apply  an  explicit  multi-grid  technique.  Both  the  unsteady  approaches  are  based  on  a dual- 
time stepping  extension  of  the  basic  time  marching  procedure  developed  for  steady  state 
computations.  Giles,  (1990)  describes  a similar  numerical  scheme  with  the  application  of  a different 
approach  for  a generic  not  even  stator-rotor  blade  count  ratio.  A recent  work  of  Dorney  et  al. 
(1997),  compares  accuracy  and  efficiency  of  different  predicting  models  with  increasing 
complexity,  ranging  from  fully  unsteady  simulations  to  loosely  coupled  methods,  for  stage 
performances  prediction.  Other  computational  studies  for  turbine  row  interaction  have  been  carried 
out  by  several  authors  among  which  Haa  et  al.  (1993),  Gallus  et  al.  (1994),  Dawes  (1994),  Fie 
(1996)  and  Walraevens  et  al.  (1998).  More  recently  von  FIoyningen-FIuene  and  Jung,  (2000)  have 
reported  the  application  and  comparison  of  different  unsteady  approaches  for  structured  grids.  The 
comparison  considers  an  explicit  method,  a time  consistent  multi-grid  and  a dual-time  step  with 
multi-grid  on  structured  grids. 

A common  drawback  of  all  the  structured  approaches  is  the  crude  way  required  to  improve 
accuracy  in  confined  region  of  rapid  gradient  which  also  implies  a costly  grid  refinement  in  a large, 
and  not  necessary,  portion  of  the  flow  domain.  A more  sensible  approach  would  instead  refine  the 
mesh  locally  only  in  regions  where  a sharp  variation  of  the  solution  is  expected.  The  first  examples 
of  unstructured  codes  for  compressible  simulations  come  from  the  external  aerodynamics  while 
more  recently  viscous  extensions  have  been  reported  by  Kwon  and  Hah,  (1995)  and  Mavriplis, 
(1995).  The  more  rational  mesh  refinement  and  the  higher  geometrical  flexibility  are  the  most 
attractive  aspects  of  the  unstructured  approaches,  allowing  complex  configurations  to  be 
represented  and  easily  handled  by  the  solution  algorithm. 

The  main  drawback  is,  undeniably,  the  high  CPU  and  memory  demand,  which  can  rapidly 
waste  the  advantages  offered  by  the  increased  geometrical  flexibility.  In  this  regard  very  few 
general  applications  have  been  documented  in  literature  on  the  use  of  unstructured  solvers  for 
unsteady  stator  rotor  interaction.  An  example  is  reported  by  Sayma  et  al.  (2000),  using  semi- 
structured  grids  with  a dual  time  step  approach.  A different  strategy  is  proposed  by  Rai  (1989) 
aiming  to  improve  the  geometrical  flexibility  through  the  use  of  patched  structured  grids. 

The  present  work  considers  the  application  of  a fully  unstructured  hybrid  solver  (HybFlow, 
Adami  et  al.  1998)  for  internal  viscous  flows.  The  multiblock  version  of  the  solver  developed  for 
turbine  rows  (Adami  et  al.,  2000)  is  considered,  because  of  the  highly  improved  performance  as 
compared  to  the  single  domain  version  of  the  code.  Moreover,  the  high  numerical  costs  involved  in 
3D  unsteady  computations  required  the  development  of  a new  parallel  single  program  multiple-data 
version  of  the  numerical  solver.  This  improved  version  of  HybFlow  is  applied  to  the  simulation  of 
the  BRITE  HP  turbine  stage  experimentally  tested  in  the  compression  tube  facility  CT3  of  the  Von 
Karman  Institute  (Denos  et  al.  1999,  2000). 
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2.  THE  GOVERNING  EQUATIONS 

The  basic  numerical  scheme  of  HybFlow  code  solves  the  compressible  Navier-Stokes 
equations  cast  in  strong  conservative  form  (Adami,  1998).  For  the  simulation  of  the  flow  field 
conditions  inside  turbine  rotor  rows,  the  conventional  condensed  formalism  has  been  adopted: 

dQ  df  dg  dh 

— + — + — + — = S (1) 

dt  dx  dy  dz 

The  governing  equations  have  been  extended  to  cope  with  a moving  frame  of  reference 
rotating  together  with  the  rotor  row.  The  same  formalism  of  governing  equations  can  be  retained  for 
both  the  fixed  and  the  rotating  frame  adding  the  centrifugal  and  Coriolis  terms  in  S , provided  that, 
when  applied  to  the  rotor  row,  the  velocity  vector  and  all  total  quantities  are  referred  to  the  relative 
frame  (Belardini  et.  al  2000).  The  perfect  gas  state  equation  is  finally  used  for  the  closure  of  the 
mathematical  model. 


3.  THE  NUMERICAL  METHOD 


The  spatial  discretization 

The  solver  HybFlow  is  specialised  for  the  computation  of  internal 
compressible/incompressible  flows  with  and  without  chemical  reactions.  A brief  description  of  the 
basic  numerical  scheme  follows,  while  more  details  can  be  found  in  Adami  et  al.  1998,  Adami  et  al. 
2000,  Belardini  et.  al.  2000.  The  spatial  discretization  is  based  on  a finite  volume  approach  for 
hybrid  unstructured  grids.  Roe's  approximate  method  is  used  for  the  upwind  scheme.  A linear 
reconstruction  of  the  solution  inside  the  elements  provides  a second  order  discretisation,  and 
monotonicity  is  ensured  through  the  TVD  slope  limiter  proposed  by  Barth,  (1991). 

Two  possible  time  accurate  discretization  are  available: 

0 The  explicit  approach;  the  time  accurate  solution  is  computed  using  a five-step  multi- 
stage Runge-Kutta  scheme,  with  the  classical  scheme. 

0 The  implicit  dual-time  stepping  discretization; 

The  dual-time  stepping  approach  adds  an  extra  numerical  time  derivative  to  the  physical 
unsteady  equation  (1): 


dQ 

dx 


^-  + R(Q)-S(Q) 
dt 


= 0 


(2) 


The  physical  time  derivative  can  be  expressed  using  a second  order  back-ward  finite 
difference: 


A Q 

At 


(*) 


phys 


3Q(k)-4Qn  +Q"-1 
2A t 


A classical  time  marching  approach  is  then  recovered  to  drive  to  convergence  the  numerical 
unsteady  term.  The  time  marching  procedure  is  based  on  the  approximate  implicit  Newton  method 
for  systems  of  non-linear  equations.  The  physical  time  derivative  is  discretized  with  second  order 
backward  approximation  and,  collecting  together  the  terms  of  the  implicit  Jacobian  matrix,  the 
classical  formalism  of  the  implicit  method  is  recovered  as  follows: 

{^+ = 
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where 


f(Q(k))= 


d 

dQ 


A Q 


At 


(k) 

phys 


+ R(Q{k))-S(Qik}) 


R 


)=- 

A Q 

(k) 

+ R(Qik))- s(Qik)) 

At 

phys 

(3) 


When  the  selected  convergence  criterion  is  satisfied  AQ'k'1  ~0  and  the  physical  solution  at  time 
level  n+1  is  updated:  Qn+ 1 =Q^ 


Stability  of  the  numerical  algorithm  is  provided  by  the  time-marching  relaxation  term 
appearing  in  the  implicit  operator  and  resulting  from  the  numerical  time  derivative  I /Ax.  The 
matrix  of  the  implicit  method  is  computed  numerically  by  finite  differences  expressions  of  the 
residual  vector  R with  respect  to  the  solution  (see  Adami,  1998),  while  the  contribution  from  the 
source  terms  is  computed  through  analytical  derivatives  of  S.  The  linear  system  stemming  from 
equation  (3)  is  solved  at  each  integration  step  by  the  same  iterative  method  GMRES  (Saad,  1994) 
used  for  the  steady  state  simulations.  The  efficient  convergence  of  the  linear  solution  requires 
preconditioning.  The  preconditioning  matrix  is  computed  performing  an  incomplete  ILU(O) 
factorisation  of  the  implicit  matrix  (Saad,  1994).  It  is  worthwhile  remembering  that  the  whole 
procedure  GMRES-ILU(O)  makes  use  of  a condensed  storage  format  which  accounts  for  matrices 
with  non-zero  elements  only. 

As  proved  in  Fig.  1-a  the  implicit  dual  time  stepping  strategy  shows  an  effective  behavior  in 
terms  of  robustness  and  stability.  In  fact  the  numerical  time  marching  converges  to  the  unsteady 
physical  solution  within  10  numerical  sub-iterations  with  about  3 orders  reduction  of  the  initial 
residuals. 


4.  THE  PARALLEL  MULTI-BLOCK  APPROACH  LOR  STAGE 
COMPUTATIONS 


The  multi-block  and  parallel  strategy 

Although  basic  solver  provides  an  accurate  description  of  internal  flows  at  different  Mach 
numbers,  when  dealing  with  complex  geometry  the  simulation  can  require  high  CPU  time  and 
memory  storage.  A significant  memory  saving  and  reduction  of  the  overall  CPU  cost  has  been 
obtained  by  means  of  different  numerical  tools.  Firstly  a multi-block  strategy  implementation 
allowed  a significant  reduction  of  memory  usage,  thanks  to  the  reduced  dimensions  of  the  linear 
system  matrix  during  the  implicit  marching  of  the  solution.  The  multiblock  approach  had  a positive 
effect  on  computing  time  also.  In  spite  of  the  loss  in  robustness  of  the  solver  due  to  the  additional 
new  internal  explicit  boundaries  created,  the  overall  computational  cost  is  considerably  reduced. 

Besides  the  domain  decomposition  technique  represented  an  essential  step,  in  view  of  the 
implementation  of  the  code  on  parallel  architectures.  The  parallel  version  of  the  code  is  based  on 
the  standard  MPI  message  passing  libraries  to  ensure  high  portability.  The  neighbouring  elements 
residing  on  the  local  memory  of  different  CPUs  require  the  explicit  activation  of  a communication 
procedure  among  processors  to  satisfy  the  physical  flow  continuity.  The  increase  in  the  number  of 
processors  leads  to  a considerable  decrease  of  computational  costs.  Nevertheless,  communications 
amongst  different  processors  represent  a computational  overhead,  which  can  grow  rapidly  with  the 
number  of  processors  involved  if  the  computational  activity  is  incorrectly  shared. 

The  present  parallel  approach  distributes  individual  portions  of  the  overall  grid  among  different 
CPUs  (data  partitioning).  In  this  case  all  processors  perform  the  same  set  of  operations  solving  the 
flow  field  inside  the  sub-domains  assigned  to  each  CPU:  the  distribution  strategy,  aimed  to 
guarantee  almost  the  same  number  of  elements  inside  each  block  and  a quite  uniform  distribution  of 
blocks  to  different  Pes,  thus  allowing  a global  load  balancing. 
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In  this  regard  the  histogram  of  Figure  1 reports  the  early  time  profiling  of  the  parallel  solver 
(percentage  of  the  different  computational  activities  with  respect  to  the  total  CPU  time),  during  the 
solution  of  a test  case.  The  test  run  on  16  processors  on  CINECA  Origin  3800,  a shared  memory 
multiprocessor  system  based  on  Processing  Elements  (PEs)  R 12000  running  at  400  MHz.  Although 
physically  distributed,  the  memory  can  be  logically  shared  by  the  hardware  able  to  address  a unified 
global  space. 

In  the  porting  of  scalar  solvers  to  parallel  architectures  a crucial  feature  in  the  estimation  of  code 
performance  is  the  communication  cost  which  is  strictly  related  to  the  structure  of  the  code  and  the 
algorithm.  Usually  for  an  implicit  solver  the  system  matrix  storage  and  inversion  require  around 
75-^80%  of  the  total  computational  time.  During  this  time  the  algorithm  does  not  require  inter- 
processor communications. 

Early  analysis  of  the  scalar  solver  performance  showed  that,  despite  the  use  of  efficient 
numerical  strategies  and  of  a compact  storage  method  (Adami  et  al.,  1998),  the  implicit  system 
required  about  88%  of  the  code  memory  demand.  Overall,  the  time  spent  for  residuals  assembly  and 
system  matrix  inversion  was  more  than  90%  of  the  total  elapsed  time.  A more  recent  profiling  of  the 
parallel  solver  shows  that  matrix  inversion  is  still  the  most  time  consuming  activity  requiring  about 
65%  of  the  global  CPU  time,  followed  by  communications  (15%),  flux  assembly  and  I/O 
management.  The  distribution  of  these  percentages  is  obviously  dependent  from  the  particular 
processor  but  the  maximum  difference  remained  below  3%  demonstrating  an  acceptable  load 
repartition.  In  view  of  these  features  of  the  solver  all  the  optimisation  strategies  performed  on  the 
parallel  code  were  mainly  focused  on  the  reduction  of  system  inversion  CPU  time  and 
communication  strategy  Among  these  the  use  of  the  most  advanced  compilers  options  ( such  as  the 
better  exploitation  of  pipeline  mechanism,  fast  improved  math  libraries  and  the  optimisation  of  the 
cache  use)  brought  a 16%  global  CPU  time  reduction  on  numerical  experiments  performed  in 
simpler  configurations.  Beside  the  Fortran  routines  managing  the  message  passage  MPI  standard 
where  implemented  with  the  introduction  of  Fortran  90  dynamic  allocation  to  reduce  the  dimension 
of  the  data  packages  to  be  sent  or  received  during  communications.  This  allowed  another  3%  CPU 
time  reduction  on  a local  4 CPUs  workstation.  A different  storage  of  system  matrix  lead  to  a 
limitation  of  global  primary  ‘cache  miss’  and  thus  a further  6%  reduction  of  the  global  CPU  time. 


Other  Activity  17-21%  (mcludiiigl/O  and 

MPI  Coinmmiicafioii 

Flux  Assembly  9-11% 


System  Matrix  Inversion  »5% 


| E3  pi«e  lilt  ■pnegmres  □laela^rn  naisgrorcs  ■ imtelte  nano  lilt  BE<-tj 

Figure  1 Time  Profiling  of  the  Solver  on  16  CPU 


To  face  the  greater  and  greater  dimensions  of  the  flow  problems  and  the  management  and  storage  of 
the  amount  of  the  data  involved,  attention  can  not  be  focused  only  on  the  optimisation  of  the  flow 
solver  itself  but  should  take  into  account  even  all  the  routines  of  the  so  called  pre  and  post- 
processing phases  which  had  to  be  revised  and  improved.  In  this  respect  Figure  2 reports  the  main 
phases  of  a multi  stage  CFD  computation  which  starts  from  the  solid  modelling  and  grid  generation 
of  the  single  stator  and  rotor  rows  using  commercial  software  (Centaur™).  The  complete  grid 
geometry  of  the  stage  is  then  assembled  in  a unique  domain  which  is  divided  in  smaller  blocks 
trough  a fully  automatic  domain  decomposition  technique  (see  Fig.  2).  Finally  the  multiblock 
domain  is  distributed  amongst  the  different  processors  and  the  processing  of  the  data  for  unsteady 
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computation  is  performed  ( Parallel  and  unsteady  treatment  in  Fig.  2).  At  this  point  the  inlet  data  set 
is  ready  for  the  CFD  simulation  which  will  produced  the  basic  flow  features  at  various  time  steps  of 
the  simulation.  Finally  the  data  produced  will  be  processed  and  synthesised  to  produce  variables 
easy  to  be  understood  and  interpreted  to  capture  the  basic  flow  physics  of  the  problem  such  as  time 
averaged  pressure  on  the  blade  surface,  unsteady  pressure  fluctuation  in  some  interesting  points  or 
for  the  comparison  with  experimental  data,  pitch  -wise  time  averages  for  secondary  flow  description 
and  movies.  With  simpler  flow  configurations  the  pre  and  post  processing  phases  where  not  critical 
but  with  the  increasing  of  the  problems  dimension  all  these  routines  had  to  be  revised  and 
improved.  In  fact  with  increasing  grid  dimensions  the  CFD  solver  is  surprisingly  not  the  crucial 
phase  of  the  whole  approach  while  the  most  memory  demanding  step  is  the  domain  decomposition 
and  the  most  CPU  intensive  step  is  the  computing  of  unsteady  information  for  the  sliding  mesh 
nterface  (see  next  paragraph).  The  reason  for  this  behaviour  can  be  explained  remembering  that 
with  increasing  grid  dimensions  a larger  CPU’s  cluster  may  be  used  to  tackle  the  computational  cost 
increase  during  the  CFD  application  having  a fully  parallel  structure.  Conversely  the  same  does  not 
hold  for  the  pre  and  post  processing  routines  which  naturally  retain  a scalar  character  required  to 
manage  the  whole  physical  domain. 

This  problems  has  been  presently  reduced,  keeping  the  same  structure  of  the  pre  processing 
routines,  with  a proper  use  of  the  memory  space  through  the  Fortran  90  standard  which  allows  a 
flexible  dynamic  memory  allocation  and  de-allocation  of  unused  matrices.  With  regard  to  the 
unsteady  treatment  the  main  problem  was  represented  by  the  length  of  the  interpolation  searching 
techniques  which  required  more  than  48  hours  to  complete  the  final  set  up  of  the  data  for  the 
unsteady  simulation  of  the  viscous  grid  made  of  1.8M  elements.  In  this  case  a more  efficient 
searching  strategy  for  the  interpolation  technique  and  a rationale  use  of  the  I/O  management  have 
been  introduced  to  obtain  a significant  reduction  of  total  required  time  (about  ten  times). 


Figure  2:  Pre  and  Post  ProcessingTreatnient 


Boundary  Conditions  Sliding  Interface 

Inlet  conditions  are  given  at  the  stator  inlet  plane  using  the  nominal  total  pressure  and  total 
temperature  profiles  as  also  the  inlet  flow  angle.  In  the  exit  plane,  hub  static  pressure  is  imposed 
and  the  radial  equilibrium  equation  is  applied  to  evaluate  the  static  pressure  distribution  assuming 
an  axially  developed  flow  field.  Solid  surfaces  are  assumed  to  be  adiabatic  with  zero  normal  mass 
and  momentum  fluxes.  Periodicity  is  ensured  with  an  exact  point  to  point  correspondence  between 
elements  belonging  to  the  periodic  boundaries.  On  the  rotor  stator  interface  a local  interpolation  has 
been  developed  to  account  for  the  relative  rows  movement  and  for  the  change  of  reference  frame. 
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Figure  3:  Interpolation  Strategy  in  Sliding  Interface 


The  grid  close  to  the  rotor-stator  boundary  is  fully  three-dimensional  and  unstructured  to  allow 
simple  and  efficient  grid  generation  in  the  inter-stage  gap.  Thus,  elements  with  triangular  and 
quadrilateral  faces  from  both  stator  and  rotor  grids  need  to  be  interfaced  at  varying  angular 
positions  (Figure  3).  In  such  situations,  the  map  of  neighbouring  elements  and  the  definition  of  a 
simple  matching  strategy  for  the  two  rows  poses  serious  difficulties,  also  because  of  the  continuous 
displacement  of  the  two  moving  blades  during  the  unsteady  computation.  This  geometrical  problem 
has  been  overcome  by  the  implementation  of  a crude,  but  general  interpolation  scheme  between  the 
two  rows,  as  summarised  in  Figure  3.  The  interpolation  strategy  guarantees  the  continuity  of  static 
variables  (density,  pressure  and  temperature),  while  velocity  components  are  discontinuous, 
according  to  the  relationships  between  absolute  and  relative  frames.  Although  the  present  method  is 
not  strictly  conservative,  the  computations  revealed  that  the  mass  imbalance  is  always  well  below 
1%.  A substantial  reduction  in  the  computing  time  to  complete  the  full  unsteady  simulation  (made 
up  of  several  periodic  revolutions),  can  be  obtained  with  the  storage  of  all  information  pertaining  to 
the  sliding  interface  interpolation  in  a look-up  table  accompanying  the  grid  file.  In  this  way  the 
parameters  for  corresponding  time  periodic  positions  can  be  computed  once  for  all  time  steps, 
avoiding  time-consuming  recalculations. 

Grid  Sequencing  Strategy 

In  a CFD  simulation  usually  an  initial  guess  solution  is  used  to  start  the  computation.  This  rough 
solution  is  updated  till  the  solver  reaches  the  final  exact  numerical  solution  of  the  physical  problem. 
When  a complex  domain  and  different  physical  aspects  are  involved,  such  as  coolant  injection 
mixed  with  the  coupling  of  rotating  and  stationary  parts,  rarely  this  initial  guess  solution  can  be 
chosen  close  to  the  real  one  and  thus  stability  problems  may  arise  while  several  time  steps  can  be 
required  before  the  final  convergence  is  obtained. 

The  number  of  this  numerical  steps  has  been  strongly  reduced,  mainly  in  the  viscous 
computation,  using  the  grid  sequencing  strategy.  The  basic  idea  beyond  this  method  is  the  creation 
and  use  of  a family  of  grids  with  progressively  increasing  element  number.  An  unsteady  unviscous 
computation  can  be  performed  on  the  coarser  grid  to  capture  the  fundamental  flow  physics  or  the 
basic  frequency  spectrum  of  the  unsteady  pressure  fluctuations  which  are  both  almost  independent 
from  viscous  phenomena.  In  the  present  application  the  coarse  unviscid  grid  is  made  up  of  360000 
elements,  of  which  24000  on  blade  surface  and  5700  on  the  sliding  surface  dividing  the  NGV  from 
the  rotor  rows. 
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Figure  5:  Refined  Grid  for  Viscous  Computations 


The  final  converged  numerical  solution  for  unviscid  simulation  still  requires  a large  number  of 
steps  but  the  CPU  time  is  reduced  in  view  of  the  smaller  dimensions  of  the  grid.  Once  the  solution 
is  obtained  on  the  coarse  grid  a viscous  and  turbulent  computations  can  be  carried  out  on  the  second 
refined  grid  (Figure  4a)  to  capture  the  effect  boundary  layers  development  on  the  performance  and 
efficiency  of  the  stage.  The  initial  guess  solution  can  be  represented  by  the  unviscid  solution 
obtained  on  the  coarse  grid  which  can  be  transferred  to  the  refined  simulation  through  an  automatic 
grid  adaptation  technique.  A general  unstructured  3D  routine  has  been  ad  hoc  developed  to 
interpolate  the  solution  from  the  coarse  to  the  second  refined  grid  which  is  reported  in  Figure  4.  A 
further  refinement  step  has  then  performed  using  a third  viscous  grid  shown  in  Figure  5.  Grid 
refinement  is  evident  in  the  radial  section  at  midspan  (Figure  5)  where  10  prismatic  layers  have 
been  introduced  all  around  the  NGV  and  rotor  blade  surfaces  to  capture  the  basic  features  of  the 
boundary  layer.  Node  elements  have  been  also  added  in  the  rotor  vane  passage  flow  and  mainly  in 
the  coolant  injection  area  (slide  (a)  in  figure  5)  of  the  NGV  which  is  thought  to  be  responsible, 
together  with  the  shock  detaching  from  the  TE  of  the  stator  blade,  for  the  unsteady  pressure 
fluctuations  measured  on  the  rotor  blade.  A strong  grid  refinement  is  also  applied  all  over  the  pitch- 
wise  direction  ahead  the  rotor  blade  (slide  (b))  in  order  to  accurately  capture  the  effect  of  the  NGV 
wake  for  all  relative  positions  of  the  two  rows.  Numerical  computations  performed  demonstrated 
that  this  refinement  requires  a soft  smoothing  toward  the  trailing  edge  of  the  rotor  to  avoid  the 
effect  of  the  grid  sensitivity  of  the  solution  on  the  different  sizes  of  neighbouring  cells  which  can 
strongly  affect  the  accuracy/stability  of  the  flow  solution.  Despite  the  refinement  the  global  grid  has 
a still  feasible  number  of  elements  (1.8  Ml  elements,  41000  on  solid  surfaces  and  20000  on  the 
sliding  plane)  owing  to  the  high  flexibility  of  the  unstructured  approach. 


5.  APPLICATION  TO  THE  3D  VKI  ANNULAR  STAGE 


Stage  description 

The  present  investigation  is  carried  out  on  the  first  stage  of  the  HP  turbine  with  prismatic  stators 
and  fully  three-dimensional  rotor  selected  in  the  frame  of  the  TATEF  project.  The  project  is  devoted 
to  gain  a better  understanding  of  the  stator-rotor  interaction  phenomena  for  transonic  annular 
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cascades.  A brief  description  of  the  cascade  is  presented  here,  while  more  details  about  the 
experimental  test  rig  and  the  data  can  be  found  in  Denos  et  al.  1999. 

A detailed  measurement  campaign  provides  data  in  several  different  operating  conditions.  In  the 
present  work  nominal  unsteady  flow  field  conditions  will  be  analysed.  The  main  mid-span 
parameters  of  the  test  case  at  the  nominal  pressure  ratio  and  for  the  nominal  reference  Reynolds 
Number  (~  106 ) are: 

Mis(exit  NGV)=1 .08,  M(exit  rotor)=0.42. 

Cooling  air  injection  is  provided  in  the  experimental  configuration  at  the  nozzle  trailing  edge 
through  a pressure  side  cut  of  the  blade.. 

Numerical  simulation 

The  3D  unstructured  grid  topology  used  to  perform  the  in  viscid  numerical  simulation  is  shown 
in  figure  4.  Applying  a slight  stretching  (about  0.8%)  to  the  rotor  vane  geometry,  an  exact  pitch 
ratio  is  obtained  using  2 stator  blades  and  3 rotor  blades.  This  modification  to  the  geometry  was 
positively  tested  by  Michelassi  et  al.  (1999). 

The  unviscid  and  viscous  (first  refinement)  computations  have  been  performed  using  the 
parallel  version  of  the  code  on  a cluster  of  four  DEC  ALPHA-XP 1000-666  Mhz  workstations.  The 
coarse  unviscid  grid  has  been  divided  using  24  blocks  equally  distributed  amongst  the  4 CPUs  with 
a memory  request  of  the  implicit  procedure  of  55  Mb. 

The  unsteady  periodic  solution  has  been  reached  starting  from  the  rest  condition  and  keeping  the 
rotor  blades  fixed.  The  starting  static  pressure  distribution  has  been  established  using  a linear  law 
between  the  inlet  prescribed  total  pressure  and  the  outlet  measured  static  pressure.  During  the 
computation  the  rotating  speed  of  the  blades  is  gradually  increased  until  the  nominal  velocity  of 
6500  RPM  is  reached  in  a physical  period  (three  rotor  passages).  The  physical  time  step  is 
computed  by  dividing  the  natural  periodic  angle  in  a suitable  number  of  steps  able  to  guarantee  both 
stability  of  the  numerical  procedure  and  an  acceptable  time  resolution  of  the  unsteady  frequencies. 
The  physical  time,  required  to  complete  an  entire  period  of  the  whole  domain,  can  be  easily 
computed  from  the  ratio  a/(0=  2.88  * HR1  .v  where  (0  is  the  rotational  speed  and  a is  the  periodic 
angle  (16.744  Deg).  The  nominal  total  quantities  are  used  as  boundary  conditions  at  the  NGV  inlet 
assuming  flat  profiles  for  both  TO  and  P0.  For  the  coolant  jet,  inlet  conditions  are  set  in  order  to 
match  a global  mass  flow  rate  as  observed  in  the  experimental  tests.  At  the  outlet  of  the  stage  the 
integral  value  of  the  nominal  experimental  outlet  static  pressure  is  imposed. 

In  this  regard  one  hundred  steps  per  subdivision  for  the  periodic  pitch  seemed  enough  to 
guarantee  the  accuracy  of  the  unsteady  solution  and  to  capture  the  first  main  frequencies  of  the 
unsteady  flow.  The  numerical  time  marching  convergence  of  each  unsteady  physical  solution  can  be 
achieved  within  10  numerical  sub-iterations  with  a reduction  of  the  initial  residuals  larger  than  4 
orders  of  magnitude.  The  unsteady  flow  field  converges  to  a pattern  periodic  in  time  after  five 
complete  passages  of  2 stators  and  3 rotors  vanes.  The  viscous  solution  has  been  obtained  starting 
from  the  solution  obtained  from  the  in  viscid  computation  using  the  grid  sequencing  strategies  and 
the  ad  hoc  interpolation  technique  previously  mentioned. 

The  computations  have  been  performed  with  the  same  solver  but  using  16  processors  on 
CINECA  Origin  3800.  The  grid  has  been  divided  in  about  400  blocks  distributed  amongst  the  16 
processing  elements  and  two  hundred  steps  have  been  chosen  to  complete  the  periodic  passage.  The 
number  of  sub  iteration  is  10  while  the  unsteady  flow  field  converged  to  a new  periodic  pattern  in 
time  after  3 complete  passages. 

Results  and  discussion 

Inviscid  results 

Once  the  periodic  solution  is  achieved  for  the  unsteady  numerical  computation,  an  overall  mass 
error  less  than  1%  has  been  observed  from  the  stage  inlet  and  outlet  sections. 

Hub  and  tip  static  pressure  profiles  are  measured  in  the  test  rig  at  a station  placed  0.035  axial 
chords  downstream  the  nozzle  trailing  edge,  using  both  Kulite  transducers  and  static  pressure  taps 
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(Denos  et  al.,  1999).  The  instantaneous  numerical  pressure  profiles  have  been  time  averaged  both  at 
hub  and  tip  and  compared  against  measurements  in  Figure  6 for  a NGV  blade  pitch. 


a)  hub  wall 


b)  tip  wall 


Figure  6:  hub  and  tip  time-averaged  pressure  profiles 


a)  Averaged  rotor  pressure  - 50%  span  b)  Averaged  rotor  pressure  - 15%  span 


c)  Averaged  rotor  pressure  - 85%  span 


Figure  7:  Time  averaged  Pressure  Profiles 


Moving  from  0 to  1 along  the  abscissa,  the  measuring  points  move  from  the  pressure  side 
towards  the  suction  side  of  the  nozzle  blade.  Fig.  6 indicates  a quite  correct  estimate  of  the  average 
inter-stage  working  pressure  in  the  numerical  simulation.  The  steep  pressure  rise  observed  for  the 
phase  value  around  0.2  can  be  explained  by  the  over  estimation  of  vane  trailing  shock  intensity  in 
the  numerical  computation  caused  by  the  absence  in  the  calculation  of  viscous  effects.  This 
behaviour  is  more  pronounced  near  the  hub  surface  where  the  inviscid  assumption  is  expected  to 
give  a less  accurate  prediction  in  view  of  the  more  relevant  role  played  by  the  boundary  layer 
thickness  on  the  lower  wall  contouring.  The  time-averaged  pressure  profiles  on  the  rotor  blade  are 
reported  in  Figure  7,  for  three  radial  sections  at  50%,  15%  and  85%  of  the  blade  span.  In  these 
figures  static  pressure  over  stage  inlet  total  pressure  is  plotted  against  curvilinear  abscissa  on  the 
rotor  surface  (positive/negative  values  indicate  respectively  locations  on  the  suction/pressure  side). 

The  agreement  with  experiments  is  satisfactory  on  the  pressure  side,  while  on  the  suction  side 
and  close  to  s/smax=0.2,  the  flow  acceleration  is  slightly  under  predicted.  This  might  indicate  either 
a three-dimensional  effect  or  a slight  underestimation  of  the  incidence  (see  Michelassi  et  al.  1999, 
2001).  The  pressure  rise  and  the  shock  location  is  correctly  captured  for  s/smax=0.6. 

At  the  15%  section  close  to  the  hub  the  computed  pressure  profile  agrees  closely  with  the 
experiments  except  for  s/smax=0. 1 on  the  suction  side.  In  this  location  it  is  argued  that  the  strong 
shock  detaching  from  the  stator  trailing  edge  interacts  with  the  rotor  blades  suction  side  giving  the 
observed  pressure  rise  on  the  rotor  surface. 

The  instantaneous  views  of  wall  pressure  patterns,  indicate  that  the  flow  on  the  LE  of  the 
suction  side  may  periodically  undergo  a strong  acceleration,  followed  by  a recompression  shock  in 
correspondence  of  some  relative  positions.  For  other  phase  angles,  the  flow  pattern  is  quite 
different,  and  the  same  acceleration/compression  pattern  is  not  detected.  A similar  pattern  can  be 
detected  close  to  tip,  although  slightly  smeared  with  respect  to  the  15%  section.  This  effect  is 
responsible  for  the  steep  pressure  rise  shown  by  the  hub  averaged  pressure  comparison  of  Figure  7. 

The  near  tip  region  pressure  field  of  figure  7 reveals  a stronger  acceleration  of  the  flow  on  the 
suction  side  close  to  the  rotor  trailing  edge  (s/smax=0.6).  This  over  expansion  is  followed, 
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downstream  the  minimum  value,  by  the  pressure  rise  typical  of  a recompression  shock,  which  is 
more  pronounced  than  in  the  experiments. 

The  absence  of  smoothing  effect  of  the  fluid  viscosity  on  the  pressure  rise  in  a shock  wave 
seems  to  be  partially  responsible  for  the  overestimation  of  pressure  unsteadiness.  Moreover,  the 
flow  pattern  at  85%  span,  close  to  the  trailing  edge  is  affected  by  the  tip  clearance,  as  indicated  by 
the  stage  simulations  by  Michelassi  et  al.  (2001).  As  a consequence  a significant  mass  can  flow 
from  pressure  to  suction  side,  thereby  decreasing  the  blade  load  and,  ultimately,  increasing  the 
minimum  pressure  level  on  the  suction  side  TE. 

Figure  8 reports  the  pitch-wise  unsteady  averaged  quantities  obtained  with  the  inviscid 
simulations.  The  axial  flow  angle  distribution  at  the  rotor  blade  exit  is  compared  with  both 
experimental  data  and  the  results  obtained  from  a NISRE  computation  (from  Denos  et  al.  1999)  in 
which  losses  are  evenly  distributed  along  the  rotor  span.  The  agreement  between  the  NISRE 
approach  and  the  unsteady  averaged  profdes  is  good.  This  is  not  surprising  since  both  approaches 
do  not  predict  any  sort  of  secondary  flow.  Actually  the  waving  shape  of  the  measured  exit  flow 
angle  is  due  to  the  combined  effect  of  various  secondary  effects,  mainly  the  presence  of  the  blade 
passage  vortex  and  tip-hub  horse  shoe  vortices. 

Figure  8a  also  shows  that  the  secondary  flows  penetrate  quite  close  to  the  rotor  mid-span,  as 
proved  by  the  deviation  from  the  nominal  exit  flow  angle  (12-deg  at  mid-span).  Figure  8b  compares 
the  experimental  time  and  pitch-wise-averaged  total  temperature  distribution  with  the  unsteady 
numerical  results.  For  radial  position  above  0.4  the  total  temperature  profile  is  followed  closely,  but 
the  absence  of  secondary  effects  in  the  calculations  provokes  a nearly  constant  total  temperaure 
distribution  which  is  far  from  experiments  close  to  the  root  portion  of  the  blade.  On  account  of  what 
seen  for  the  exit  flow  angle,  this  disagreement  is  not  surprising.  Still  the  unsteady  simulations 
demonstrate  that  the  mid-span  value  of  the  total  pressure  is  remarkably  well  reproduced,  thereby 
proving  the  existence  of  a nearly  two-dimensional  flow  portion  around  the  rotor  blade.  The 
agreement  with  the  NISRE  calculations  are  good,  although  both  are  not  able  to  reproduce  the 
correct  total  temperature  profile  which,  as  seen  for  the  exit  flow  angle,  is  strongly  affected  by 
secondary  flows  which,  while  governing  the  exit  flow  angle,  can  alter  the  local  extraction  of  energy. 
Globally  the  integral  value  of  the  inviscid  total  temperatures  is  in  good  agreement  with  the 
experimental  data  and  NISRE  estimate,  allowing  a reasonable  forecast  of  energy  exchanges.  Near 
hub,  where  secondary  flows  and  wall  boundary  layers  are  well  developed,  the  accuracy  of  the 
unsteady  computation  is  similar  to  that  obtained  with  the  NISRE  approach. 


a)  Averaged  Alfa  Flow  angle 


b)  Averaged  Total  Temperature  c)  Averaged  Total  Pressure 

Figure  8:  Pitch-wise  Averaged  Quantities 


The  total  pressure  profiles  shown  in  figure  8c  allow  the  same  conclusions  to  be  drawn.  The  time 
averaged  values  stemming  from  the  unsteady  inviscid  calculations  follow  the  experimental  results 
with  better  accuracy,  especially  in  the  lower  part  of  the  blade  in  spite  of  some  inaccuracies 
connected  to  the  inviscid  approximation.  The  capability  of  the  unsteady  approach  to  model  the  time 
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dependent  flow  should  be  checked  by  the  analysis  of  the  frequency  spectrum.  Instantaneous 
pressure  has  been  monitored  for  several  representative  points  located  along  rotor  and  stator  blade 
surface  (Denos  et.  al.,  1999,  Figure  18-a). 


a)  Rotor  Gauge  Position 


Frequency  (Hz) 

b)  Rotor  LE  gauge  4 


c)  Rotor  TE  gauge  12 


Figure  9:  Frequency  Spectrum 


The  frequency  analysis  for  two  points  located  near  LE  and  TE  of  rotor  blade  on  suction  surface 
are  reported  in  Figure  9.  For  both  gauges  the  main  harmonic  observed  is  the  stator  blades  passing 
frequency.  These  basic  frequencies  of  the  spectrum  (figures  18-b  and  18-c)  can  be  computed  from 
the  speed  rotation  (6500  RPM)  and  the  blade  count  ratio  (2:3)  leading  to  a main  peaks  at/=4.66  kFIz 
and  the  corresponding  multiples  (9.32  kHz,  13.98  kHz  ...).  The  remarkable  amplitude  of  the  first 
harmonic  at  LE  is  due  both  to  the  not  uniform  exit  flow  field  from  the  stator  blades  (mainly  due  to 
the  presence  of  vane  TE  shock  impinging  on  the  rotor)  and  to  the  small  gap  existing  between  the 
stator  TE  and  the  rotor  LE.  The  amplitude  of  the  pressure  discontinuity  varies  considerably 
depending  on  the  location  on  the  rotor  blade.  The  maximum  amplitude  is  at  gauge  4 because 
directly  heat  by  the  shock  wave.  Flowing  toward  the  TE,  the  pressure  fluctuations  are  strongly 
reduced  (two  orders)  by  the  action  of  the  vane  acceleration  effect,  as  can  be  observed  in  figure  9-c. 


Figure  10:  Unsteady  Pressure  Fluctuations 


Figure  10  shows  the  unsteady  pressure  fluctuations  Ap/P0\  = {p- pavemged)/P0\  for  two  rotor 

blade  locations  at  mid-span.  The  two  selected  points  lay  on  the  suction  side  at  approximately  20% 
axial  chord  (again  point  4 in  Figure  9)  and  on  the  leading  edge  (point  14).  On  the  axis  the  phase 
spans  from  0 to  2 while  a complete  periodic  passage  is  completed  (e.g.  2 NGV  vanes  and  3 rotor 
blades).  The  experimental  data  have  been  described  by  Denos  et  al.,  2000,  who  proposed  a detailed 
interpretation  of  the  flow  in  the  stator/rotor  gap.  In  the  experiments  pressure  traces  on  gauge  4 
clearly  indicate  high  amplitude  fluctuations  in  the  rotor  LE  region,  related  to  the  passage  of  the  vane 
TE  fish  tale  shock  (periodically  oscillating  around  the  axial  direction  from  about  +20°  to  -20°).  In 
the  same  Figure  two  time-accurate  computations,  are  compared  with  the  experimental  unsteady 
profile,  obtained  doubling  the  number  of  time  steps  during  the  unsteady  simulation  and  the  accuracy 
seems  to  improve  with  the  increasing  time  resolution.  The  basic  feature  of  the  instantaneous  signal 
are  reproduced  by  the  numerical  results  even  if  peak  amplitudes  of  pressure  fluctuations  are  under 
predicted  mainly  at  gauge  4.  Anyway  in  this  point  the  under  estimation  behaves  coherently  with  the 
pressure  disagreement  already  observed  at  the  same  station  in  the  averaged  profiles  of  Figures  7. 
Moreover,  in  this  location  the  stator  wake  hits  the  rotor  blade  almost  tangentially  (Michelassi  et  al., 
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1999).  This  could  provoke  a complex  interaction  between  the  rotor  boundary  layer  and  the  stator 
wake,  which  is  clearly  governed  by  viscous  effects. 

Viscous  results 

The  time-averaged  static  pressure  profiles  on  the  rotor  blade  are  reported  in  Figure  1 1 against 
curvilinear  abscissa,  for  the  same  radial  sections  of  Figure  7.  The  agreement  with  experiments  at 
50%  span  on  pressure  side  confirms  the  in-viscid  results.  On  the  suction  side  the  under  prediction  of 
flow  acceleration  is  still  present  confirming  the  slight  underestimation  of  the  incidence  but  the 
pressure  distribution  of  the  experimental  results  near  the  shock  location  is  more  accurately 
followed. 


a)  Averaged  rotor  pressure  - 50%  span 


b)  Averaged  rotor  pressure  - 15%  span  c)  Averaged  rotor  pressure  - 85%  span 


Figure  11:  Time  averaged  Pressure  Profiles 


As  expected  in  the  viscous  computation  the  up-down  shape  shown  in  the  in  viscid  assumption 
disappeared  due  to  the  presence  of  smoothing  effect  of  the  fluid  viscosity  on  the  level  of 
acceleration  and  recompression  in  correspondence  of  some  relative  positions  of  the  two  rows.  Also 
the  strong  acceleration  of  the  flow  on  the  suction  side  close  to  the  rotor  trailing  edge  at  tip  and  hub 
disappeared  (Figures  11a  and  lib),  according  to  the  more  accurate  physical  assumption  and  the 
agreement  with  experimental  results  is  improved. 

Phenomena  such  secondary  flow  development  inside  the  stator  blade  and  their  evolution  in  the 
rotor  blade,  depending  on  the  relative  blade  position  between  rotor  and  stator  and  the  blade  wakes 
development  have  a deep  impact  on  loss  profile  and  distribution  and  thus  on  the  whole 
turbomachine  performances. 

The  viscous  approximation  can  provide  an  useful  tool  for  the  basic  understanding  of  such 
phenomena.  In  this  regard  in  Figure  1 1 .b  the  entropy  contours  at  the  stator  exit  are  shown  in  axial 
plane  1,  situated  at  10  % of  the  axial  chord  downstream  the  stator  TE.  The  stator  wake  is  well 
outlined  by  the  high  loss  region  in  the  middle  of  the  area. 

As  can  be  noticed  the  mid-span  flow  is  essentially  2D,  most  of  the  loss  being  confined  in  the 
blade  wake.  Near  hub  and  tip  other  high  loss  areas  (1)  can  be  underlined  due  to  the  development  of 
wall  boundary  layers  in  the  stator  channel.  The  extension  and  thickness  is  small  in  comparison  to 
the  radial  and  pitchwise  extension  of  the  blade,  limiting  their  influence  in  the  regions  very  close  to 
solid  surfaces.  As  expected  the  boundary  layer  growth  doesn’t  seem  to  affect  the  core  of  the  flow, 
which  remains  almost  undisturbed. 
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d)  Entropy  Contours  - Plane  3 


Figure  12:  Entropy  Contours 

In  Figure  12  e the  entropy  contours  at  50%  span  are  shown  and  the  TE  wakes  can  be  clearly 
identified.  These  contours  confirm  that  the  entropy  production  in  the  NGV  is  essentially  provoked 
by  the  SS  boundary  layer  that  grows  until  the  final  part  of  the  blade,  detaches  in  correspondence  of 
the  stator  TE  region,  diffuses  and  finally  impinges  on  the  rotor  blade.  On  the  pressure  side  the 
boundary  layer  is  almost  invisible  till  the  LE  is  reached,  due  to  favourable  pressure  gradient 
experienced  by  the  flow  field  in  most  of  wall  surface.  On  this  side  of  the  blade  the  more  relevant 
loss  seems  to  be  represented  by  the  mixing  process  between  coolant  injection  and  the  main  flow 
downstream  the  TE,  responsible  for  the  second  high  loss  stripe  of  the  wake.  The  shape  of  the  wake 
from  the  NGV  is  almost  independent  from  the  relative  position  between  the  two  rows  while  losses 
in  the  main  core  of  the  flow  are  not  evident.  When  the  wake  moves  toward  the  rotor  blade  the  effect 
on  the  flow  can  be  different  depending  on  the  relative  position  of  the  two  rows.  For  instance  in 
figure  12.e  the  wake  impinges  on  the  LE  of  the  blade  and  is  split  in  two  different  branches,  the  first 
one  on  the  PS  and  the  second  on  the  SS.  After  a half  period  the  wake  impinges  directly  on  the  blade 
PS  effecting  primarily  the  SS  flow  on  the  LE. 

In  figure  12.d  the  entropy  carpet  plot  is  shown  for  reference  plane  3,  located  downstream  the 
rotor  blade  TE  for  tangential  position  corresponding  to  that  of  fig  12.e.  In  this  plane  high  loss 
regions  can  be  clearly  identified  corresponding  to  the  rotor  blade  wake.  Near  the  tip  entropy 
production  bubbles  can  be  detected  according  to  the  combined  effect  of  the  tip  leakage  vortex  and 
shroud  boundary  layer  that  are  both  contributing  to  the  secondary  flows  in  this  region. 

The  ad  hoc  grid  refinement  performed  for  viscous  computations  in  this  area  were  necessary  for 
the  more  accurate  description  of  the  viscous  boundary  layer  profile  and  the  velocity  distribution. 
The  shape  and  the  distribution  of  the  high  loss  region  appears  to  be  quite  dependent  on  the  relative 
row  positions  of  the  stator/rotor  rows.  Deep  discussion  of  the  results  could  be  possible  but  is  out  of 
the  scope  of  this  paper  intended  to  report  mainly  the  challenge  and  the  effort  required  in  solving  the 
unsteady  flow  field  in  real  stage  configuration. 

6.  CONCLUSIONS 

A 3D  unsteady  solver  has  been  discussed  in  the  viscous  and  in-viscid  assumption  through  the 
comparison  against  experiments.  The  peculiar  aspect  of  the  present  application,  compared  to 
documented  works,  is  given  by  the  completely  hybrid-unstructured  nature  of  the  approach.  This 


Plane  i 


b)  Entropy  Contours  - Plane  I 


e)  Entropy  Contours  - Midspan 


Plane  2 


c)  Entropy  Contours  - Plane  2 


f)  Entropy  Contours  - '/Period 
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feature  allows  an  easy  and  flexible  representation  of  the  stage  in  view  of  future  more  detailed  flow 
investigations,  as  required  by  viscous  simulations  or  shock  interactions,  and  of  more  complex 
geometry  as  required  by  cooling  systems  or  by  row  leakages. 

The  use  of  unstructured  grids  posed  new  difficulties  for  the  matching  of  the  stator/rotor  rows, 
which  have  been  here  overcome  by  a crude  but  general  and  automatic  linear  interpolation  strategy, 
but  on  the  other  hand  offers  a great  flexibility  in  the  grid  generation  and  modelling  with  substantial 
reduction  in  the  grid  points  demand.  The  potential  stator  rotor  interaction  has  been  reproduced  with 
a satisfactory  accuracy  demonstrating  that  the  unsteady  Euler  approach  allows  a more  realistic 
description  in  comparison  to  steady  state  computations  of  the  flow  pattern,  especially  in  the 
presence  of  physical  phenomena  such  unsteady  shock  interaction  which  can  not  be  accounted  by  in 
a steady  assumption. 

In  order  to  have  a deeper  insight  of  phenomena  concerning  the  periodic  evolution  of  the  NGV 
wakes,  tip  vortex  structures  and  the  effect  of  unsteadiness  on  the  global  stage  efficiency,  a more 
accurate  analysis  , based  on  full  viscous,  calculation  was  required.  Viscous  calculation  on  the  other 
hand  posed  serious  problems  for  the  management  of  the  data  file  and  the  CPU  usage.  To  accurately 
reproduce  those  viscous  effects  more  stringent  requirements  were  imposed  on  the  grid  cell 
dimension  especially  near  boundary  layers  in  the  solid  surfaces  at  hub  tip  and  blade  profile.  The 
grid  dimensions  typical  of  a realistic  viscous  computations  (about  1.8  ML  node  elements)  created 
difficulties  for  the  management  and  storage  of  the  input  and  output  data  required  or  produced 
during  the  processing  of  the  solution  and  compelled  the  complete  reconsideration  of  all  the  phases 
involved  in  the  CFD  solution  and  the  optimisation  of  all  the  pre  and  post  processing  routines.  The 
most  advanced  programming  tools  such  Fortran  90  dynamic  memory  allocation  and  de-allocation  of 
unused  matrices  were  necessary  for  the  full  storage  of  geometry,  exceeding  the  capacity  of  local 
workstation.  Efficient  interpolation  techniques  and  a rationale  use  of  the  I/O  management  were 
necessary  to  reduce  the  length  of  the  grid  sequencing  technique  and  all  the  post  processing  phases 
of  the  results  management.  In  conclusion  great  effort  have  been  successfully  carried  out  on  the 
different  aspects  of  the  calculation  to  face  correctly  and  profitable  the  challenge  to  compute  and 
investigate  unsteady  flow  in  a full  turbine  stage. 

ACKNOWLEDGMENTS 

The  present  research  was  carried  out  under  contract  for  the  European  Commission  as  part  of  the  Brite- 
EuRam  TATEF  (Trubine  Aero-Thermal  External  Flow)  project  (BRPR-CT97-0519).  The  authors  wish  to 
acknowledge  Dr.  Ing.  Paolo  Malfetti  and  CINECA  for  the  tests  on  the  Origin3k  and  for  the  active  support  in 
the  parallel  solver  porting. 

REFERENCES 

Adami  P.,  Martel  li  F.  and  Michelassi  V.,  2000,  “Three-Dimensional  Investigations  for  Axial  Turbines  by 
an  Implicit  Unstructured  Multi-block  Flow  Solver”,  ASME,  IGTI  TurboExpo  2000,  Munich. 

Adami  P.,  Michelassi,  V.,  Martelli,  F.,  1998  “Performances  of  a Newton-Krylov  scheme  against  implicit 
and  multi-grid  solvers  for  inviscid  flows”  A1AA  paper  98-2429. 

Adami,  P.,  1999  “Computations  for  internal  flows  with  a low-mach  preconditioned  Newton-Krylov 
scheme”,  3rd  European  Conference  on  Turbomachinery.  London. 

Barth,  T.J.,1991  "Numerical  Aspects  of  Computing  Viscous  High-Reynolds  Number  Flows  on 
Unstructured  Meshes",  A1AA  Paper  91-0721,  Jan. 

Belardini  E.,  Adami  P.,  Martelli  F.  “ Development  of  an  Unsteady  Parallel  Approach  for  3D  Stator-Rotor 
Interaction”  4th  European  Conference  on  Turbomachinery,  Fluid  Dynamics  and  Thermodynamics-  Firenze 
20-23  March  2001,  IMechE  Journal  of  Power  and  Enerev2001-  Vol.215.  n,A6. 

Dawes,  W.N.,  1992,  “The  Simulation  of  Three-Dimensional  Viscous  Flow  in  Turbomachinery 
Geometries  Using  a Solution- Adaptive  Unstructured  Mesh  Methodology”  Transaction  of  ASME  Vol.  114, 
July. 

Dawes,  W.N.,  1994,  "A  Numerical  Study  of  the  Interaction  of  a Transonic  Compressor  Rotor  Overtip 
Leakage  Vortex  with  the  Following  Stator  Blade  Row",  ASME  Paper  94-GT-156. 


37-16 


Denos  R..  Sieverding  C.H.,  Arts  T.,  Brouckaert  J.F.  and  Paniagua  G.,  1999,  "Experimental  Investigation 
of  the  Unsteady  Rotor  Aerodynamics  of  a Transonic  Turbine  Stage",  3rd  European  Conference  on 
Turbomachinery.  London. 

Denos,  R.,  Arts,  T.,  Paniagua  G„  Michelassi,  V.  and  Martelli,  F.,  2000,  “Investigation  of  the  unsteady 
Rotor  Aerodinamics  in  a Transonic  Turbine  Stage”  2000-GT-435. 

Dorney  D.J.  and  Sharma  O.P.,  1997,  "Evaluation  of  Flow  Field  Approximations  for  Transonic 
Compressor  Stages",  J of  Turbomeachinery,  Vol.  1 19. 

Gallus  H.E.,  Zeschky  J.  and  Hah  C.,  1994,  "Endwall  and  Unsteady  Flow  Phenomena  in  an  Axial  Turbine 
Stage"  ASME  94-GT- 143. 

Giles,  M.  B.,  1990,  “Stator/Rotor  Interaction  in  a Transonic  Turbine”  A1AA  J.  of  Propulsion,  Vol.  6, 
No.  5. 

Hah  C.,  Copenhaver  W.W.  and  Puterbaugh  S.L.,  1993,  "Unsteady  Aerodynamic  Flow  Phenomena  in  a 
Transonic  Compressor  Stage",  A1AA-93-1868. 

Haselbacher,  A.,  McGuirk,  J.J.,  Page,  G.J.,  1999,  “Finite  Volume  Discretization  Aspects  for  Viscous 
Flows  on  Mixed  Unstructured  Grids”  A1AA  J.,  Vol.  37,  No.  2. 

He  L.,  1996, "Time  Marching  Calculations  of  Unsteady  Flows,  Blade  Row  Interaction  and  Flutter”, 
VKI-LS  1996-05. 

Kwon,  O.J.  and  Hah,  C.,  1995,  “Simulation  of  Three-Dimensional  Turbulent  Flows  on  Unstructured 
Meshes”  AIAA  J.  Vol.33,  No.  6,  June. 

M.  Von  Hoyningen-Huene,  A.  R.  Jung,  2000,  “Comparison  of  Different  Acceleration  techniques  and 
methods  for  Periodic  Boundary  Treatment  in  Unsteady  Turbine  Stage  Flow  Simulations” 
J.  of  Turbomachinerv,  Vol,  122. 

Mavriplis  D.J.  1995,  “Three-Dimensional  Multigrid  Reynolds  Averaged  Navier  Stokes  Solver  for 
Unstructured  Meshes”  AIAA  J.,  Vol.33,  No.  3. 

Martelli,  F.  Adami,  P.  Belardini,  E,  “Unsteady  Rotor/Stator  Interaction  : An  Improved  Unstructured 
Approach”  46th  ASME  International  Gas  Turbine  & Aeroengine  Technical  Congress  June  4-7,  2001, 
N. Orleans,  Louisiana,  USA 

Michelassi,  V,  Martelli,  Denos,  R„  Arts,  T.,  Sieverding,  C.H.  1999  “Unsteady  Heat  Transfer  in  Stator- 
Rotor  Interaction  by  Two-Equations  Turbulence  Model”  J.  Turbomachinery, July. 

Michelassi,  V.,  Giangiacomo,  P.,  Martelli,  F.,  Denos,  R.,  Paniagua,  G.,  “Steady  three-dimensional 
simulation  of  transonic  axial  turbine  stage”,  offered  for  IGTI-2001. 

Pulliam  T.,H.,  Rogers,  S.,  Barth,  T.,  1996,  "Practical  Aspects  of  Krylov  Subspace  Iterative  Methods  in 
CFD",  Prog,  and  Challenges  in  CFD  Methods  and  Algorithms , Seville,  AGARD  CP-578. 

Rai,  M.,  M.,  1989,  “Three-Dimensional  navier-Stokes  Simulations  of  Turbine  Rotor-Stator  Interaction; 
Part  I-Methodology”  AIAA  J.  of  Propulsion,  Vol.  5,  No.  3 May-June. 

Roe  P.L.,  1986  “Characteristic  based  scheme  for  Euler  equations”  , Ann.  Rev.  Fluid  Mech.,  18. 

Saad,  Y.,  1994  "Krylov  Subspace  Thechniques,  Coniugate  Gradients,  Preconditioning  and  Sparse  Matrix 
Solvers",  CFD  VKI  LS  1994-05  VonKarman  Institute  for  Fluid  Dynamics. 

Saad,  Y.,  1994,  "Krylov  Subspace  Techniques,  Coniugate  Gradients,  Preconditioning  and  Sparse  Matrix 
Solvers",  CFD  VKI  LS  1994-05  VonKarman  Institute  for  Fluid  Dynamics. 

Sayma  A.  I.,  Vahdati  M.,  Sbardella  L.,  and  M.  Imregun,  2000,  “Modeling  of  the  Three-Dimensional 
Viscous  Compressible  Turbomachinery  Flows  Using  Unstructured  Hybrid  Grids”  A1AAJ.  Vol.  38,  No.  6. 

Sharma  O.P.,  Pickett  G.F.  and  Ni  R.H.,  1992,  "Assessment  of  Unsteady  Flows  in  Turbines",  J.  of 
Turbom.  Vol.  1 14. 

Slack,  D.,C.,  Whitaker,  D.,L.,  Walters,  R.,W.,  1994,  "Time  Integration  Algorithms  for  the  Two- 
Dimensional  Euler  Equations  on  Unstructured  Meshes",  AIAA  J.  Vol.  32,  No.  6. 

Venkatakrishnan,  V.,  1995,  "A  Perspective  on  Unstructured  Grid  flow  Solvers",  Technical  report  AIAA 
Paper  95-0667,  AIAA  33rd  Aerospace  Sciences  Conference,  Reno. 

Walraevens  R.E.,  Gallus  H.E.,  Jung  A.R.,  Mayer  F.J.  and  Stetter  H.,  1998,  "Experimental  and 
Computational  Study  of  the  Unsteady  Flow  in  a 1.5  Stage  Axial  Turbine  with  Emphasis  on  the  Secondary 
Flows  in  the  Second  Stator"  ASME  Paper  98-GT-254. 

Centaur™  -Grid  Generation  Software  Package  Distributed  by  CentaurSoft  (web  site  www.centaursoft.com) 


37-17 


Paper  #37 

Discussor’s  name  A.  Soulaimani 
Author  F.  Martelli 

Q:  1)  Since  the  precondition  and  its  factorization  take  most  of  the  CPU  time,  do  you  compute  it  at 
every  time  step? 

2)  You  can  save  a substantial  amount  of  computing  time  if  you  freeze  it  for  a few  time  steps  (2-5) 
steps).  I am  sure  of  the  result! 

A:  1)  Yes 
2)  Okay 


